#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "sparse_crs.h"


int N;				/**< 计算区域分段数. */
int n;				/**< 总自由度. */

int main(int argc,char *argv[]){
    clock_t t;
    
    FILE *input;
 
    struct PATTERN pA,pA1,pA2,pA3;
    struct MATRIX A,A1,A2,A3;


    int i,j;//列、行编号

    N = atoi(argv[1]);
    n = (N * N);

    t = clock();
    /// 准备系数矩阵模板.
    init_pattern(&pA, n, n, 5);
    /// Ui,j
    for (i = 0; i < N; i++)
	for (j = 0; j < N; j++)
	    set_pattern(&pA, j * N + i, j * N + i);
    /// Ui+1,j
    for (j = 0; j < N; j++)
	for (i = 0; i < N - 1; i++)
	    set_pattern(&pA, j * N + i, j * N + i + 1);
    /// Ui-1,j
    for (j = 0; j < N; j++)
	for (i = 1; i < N; i++)
	    set_pattern(&pA, j * N + i, j * N + i - 1);
    /// Ui1,j-1
    for (i = 0; i < N; i++)
	for (j = 1; j < N; j++)
	    set_pattern(&pA, j * N + i, j * N + i - N);
    /// Ui1,j+1
    for (i = 0; i < N; i++)
	for (j = 0; j < N - 1; j++)
	    set_pattern(&pA, j * N + i, j * N + i + N);
    compress_pattern(&pA);

    init_matrix(&A, &pA);
    /// Ui,j, 右端项
    for (i = 0; i < N; i++)
	for (j = 0; j < N; j++)
	{
	    add_matrix(&A, j * N + i, j * N + i, 4.0);
	}
    /// Ui+1,j, 右边界
    for (j = 0; j < N; j++)
    {
	for (i = 0; i < N - 1; i++)
	    add_matrix(&A, j * N + i, j * N + i + 1, -1.0);
    }
    /// Ui-1,j, 左边界
    for (j = 0; j < N; j++)
    {
	for (i = 1; i < N; i++)
	    add_matrix(&A, j * N + i, j * N + i - 1, -1.0);
    }
    /// Ui1,j-1, 下边界
    for (i = 0; i < N; i++)
    {
	for (j = 1; j < N; j++)
	    add_matrix(&A, j * N + i, j * N + i - N, -1.0);
    }
    /// Ui1,j+1, 上边界
    for (i = 0; i < N; i++)
    {
	for (j = 0; j < N - 1; j++)
	    add_matrix(&A, j * N + i, j * N + i + N, -1.0);

    }
    t = clock() - t;
    printf("The first matrix mat1 is assembled.\n");
    //printf ("The first assembling tooks me %d clicks (%f seconds).\n",(int)t,((float)t)/CLOCKS_PER_SEC);
    //i = print_matrix(&A);
    //测试转置
    printf("Test the transpose operator:\n");
    printf("The original matrix is:\n");
    //测试打印三元组
    print_triple(&A);
    init_pattern_trans(&pA3, n, n, n);//初始化列全为-1
    init_matrix_trans(&A3, &pA3);
    i = mat_trans(&A,&A3);
    printf("Complete transpose!\n");
    printf("The transpose matrix is:\n");
    print_triple(&A3);
    
    t = clock();
    /// 准备系数矩阵模板.
    init_pattern(&pA1, n, n, 5);
    /// Ui,j
    for (i = 0; i < N; i++)
	for (j = 0; j < N; j++)
	    set_pattern(&pA1, j * N + i, j * N + i);
    /// Ui+1,j
    for (j = 0; j < N; j++)
	for (i = 0; i < N - 1; i++)
	    set_pattern(&pA1, j * N + i, j * N + i + 1);
    /// Ui-1,j
    for (j = 0; j < N; j++)
	for (i = 1; i < N; i++)
	    set_pattern(&pA1, j * N + i, j * N + i - 1);
    /// Ui1,j-1
    for (i = 0; i < N; i++)
	for (j = 1; j < N; j++)
	    set_pattern(&pA1, j * N + i, j * N + i - N);
    /// Ui1,j+1
    for (i = 0; i < N; i++)
	for (j = 0; j < N - 1; j++)
	    set_pattern(&pA1, j * N + i, j * N + i + N);
    compress_pattern(&pA1);

    init_matrix(&A1, &pA1);
    /// Ui,j, 右端项
    for (i = 0; i < N; i++)
	for (j = 0; j < N; j++)
	{
	    add_matrix(&A1, j * N + i, j * N + i, 4.0);
	}
    /// Ui+1,j, 右边界
    for (j = 0; j < N; j++)
    {
	for (i = 0; i < N - 1; i++)
	    add_matrix(&A1, j * N + i, j * N + i + 1, -1.0);
    }
    /// Ui-1,j, 左边界
    for (j = 0; j < N; j++)
    {
	for (i = 1; i < N; i++)
	    add_matrix(&A1, j * N + i, j * N + i - 1, -1.0);
    }
    /// Ui1,j-1, 下边界
    for (i = 0; i < N; i++)
    {
	for (j = 1; j < N; j++)
	    add_matrix(&A1, j * N + i, j * N + i - N, -1.0);
    }
    /// Ui1,j+1, 上边界
    for (i = 0; i < N; i++)
    {
	for (j = 0; j < N - 1; j++)
	    add_matrix(&A1, j * N + i, j * N + i + N, -1.0);

    }
    t = clock() - t;
    /* printf ("The second assembling tooks me %d clicks (%f seconds).\n",(int)t,((float)t)/CLOCKS_PER_SEC); */
    printf("The second matrix mat2 is assembled.\n");

    init_pattern_source(&pA2, n, n);
    init_matrix_source(&A2, &pA2);

    t = clock();
    printf("Test the multiply.\n");
    double i1 = A.pat->row[n]*1.0/n/n;
    printf("The ratio of none zero number %d to size of matrix %d is:%lf\n",A.pat->row[n],(n*n),A.pat->row[n]*1.0/n/n);
    multiMatrix(&A,&A1,&A2);
    t = clock() - t;
    //i = print_matrix(&A2);
    
    printf ("The multiply tooks me %d clicks (%f seconds).\n",(int)t,((float)t)/CLOCKS_PER_SEC);
    //printf("The number of none zero is :%d\n",A2.pat->row[n]);
    
    //i = print_matrix(&A2);
    
    return 0;
    
}
